BIOL 3100 Final Project

Effects of maternal stress and gut microbiome on offspring social behaviors in long-bout stress

Project background

Essentially all individuals experience social stress, which may be acute or chronic. Individuals differ in their behavioral and physiological responses to stressors, but what factors cause these individual differences is less clear. Manipulations to the maternal microbiome and maternal stress can have sex-specific effects on offspring social behaviors and may help explain variation in behavioral responses to prolonged stress. My research project investigates how maternal environment affects differences in offspring social behavior during long bouts of stress in Siberian hamsters (Phodopus sungorus, figure 1).

Figure 1. Image of Siberian hamster Phodopus sungorus. Image attribution: Salzgeber, P. (2007). Phodopus sungorus. In Wikimedia Commons. https://commons.wikimedia.org/wiki/File:PhodopusSungorus_1.jpg

---

Research questions

Question: How does early development (i.e., maternal environment) affect whether offspring vary in their social behaviors during the course of a long-bout social stressor?
Hypothesis: Maternal stress during offspring development has sex-specific effects on consistency in individual behavior during a long-bout social stressor
Predictions:
1. Female offspring from Stress Only mothers will display increased aggression throughout a long-bout stressor.
2. Other offspring will not display consistent behaviors throughout the long-bout stressor.

Methods

Maternal treatment

Pregnant females (n=34) were exposed to one of four treatments for ten days during pregnancy: control, antibiotic only, stress only, antibiotic + stress (figure 2)

Figure 2. Maternal treatment graphic depicting the four treatment groups.

Offspring behavioral trials

Juvenile offspring (n=14 males, 11 females) were exposed to a 15 minute social stressor following a resident-intruder paradigm. The intruder was placed into the focal individual’s (i.e., offspring’s) cage for 15 minutes, and offspring behaviors towards the intruder were scored (figure 3).
The first 5 minutes (t0-t5) and the last five minutes (t10-t15) of behaviors were scored and analyzed.

Figure 3.Offspring in late adolescence (51-56 days old) were exposed to a resident-intruder trial as above. The intruder was placed in the focal individual’s cage for 15 minutes. Time scored is shown by red line.

Four behaviors were scored: attack, investigation, escape, and grooming (figure 4).

Figure 4. Actual screenshots of video data depicting the four social behaviors of interest.

Data cleaning

Below is the code used to clean the and prepare the data for analyses. Several columns were removed from the data frame because they were not used in this project (these were columns imported into the excel sheet from the original project in 2019).

# Load data, clean names, view structure ####
hammies <- read_csv("../final_project/IU106_Social Behavior First and Last Five Minutes_Formatted_6-6-23.csv")
hammies <- janitor::clean_names(hammies)
names(hammies)
str(hammies)

# Remove columns that were not used in this project ####
columns_removed <- c("id", "aggression_number_assigned", "comp1", "comp2", "comp3", 
                     "aggression_score_freq", "aggression_score_duration")
hammies <- dplyr::select(hammies, -columns_removed)

# Create trial timing column, set appropriate columns as factors ####
tidy_hammies <- hammies %>% 
  mutate(trial_timing = as.factor(first_or_last_five_minutes)) %>% 
  mutate(treatment = as.factor(maternal_treatment)) %>% 
  mutate(juvenile_id = as.factor(juvenile_id)) %>% 
  dplyr::select(-first_or_last_five_minutes)

# Fix spelling of all "occurence" columns ####
tidy_hammies <- tidy_hammies %>% 
  mutate(attack_number_of_occurrences = attack_number_of_occurences) %>%
  mutate(chase_number_of_occurrences = chase_number_of_occurences) %>%
  mutate(groom_number_of_occurrences = groom_number_of_occurences) %>%
  mutate(investigation_number_of_occurrences = investigation_number_of_occurences) %>%
  mutate(jump_number_of_occurrences = jump_number_of_occurences) %>% 
  dplyr::select(-contains("occurences"))

str(tidy_hammies)

Modeling Juvenile Behavior

We analyzed the number (i.e., frequency) and total duration of behavioral bouts in the first and last five minutes (t0-t5 vs. t10-t15).
For each model, we first tested for a three-way interaction between maternal treatment, offspring sex, and trial timing. If no significant three-way interaction was detected, we tested for a two-way interaction between offspring sex and maternal treatment. Only significant interactions were retained in the model; non-significant interaction terms were removed.

Summary of models

All models had two random effects: intruder identity and resident (offspring) identity, to control for individual differences in behavior. All frequency models (GLMM’s) used the bobyqa optimizer from the lme4 package for model fitting.

Table 1. Summary of models run for each response variable.

Response Variable	Distribution	Fixed Effects
Attack Frequency	Poisson, Log link function	Offspring sex, maternal treatment, trial timing, SI- CORT, offspring weight, offspring sex*maternal treatment interaction
Attack Duration	Gaussian, Identity link function	Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight
Investigation Frequency	Poisson, Log link function	Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight, offspring sexmaternal treatmenttrial timing
Investigation Duration	Gaussian, Identity link function	Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight
Escape Frequency	Poisson, Log link function	Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight, offspring sexmaternal treatmenttrial timing
Escape Duration	Gaussian, Identity link function	Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight
Grooming Frequency	Poisson, Log link function	Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight
Grooming Duration	Gaussian, Identity link function	Offspring sex, maternal treatment, trial timing, SI-CORT, offspring weight, offspring sex* maternal treatment

Agression frequency

Results

Our first model analyzes the effects of maternal treatment, offspring sex, and video timing on the frequency of aggressive behaviors. There was no significant three-way interaction, but we did detect a significant interaction between maternal treatment and offspring sex. Below is the code for model A.1.

modelA.1 <- tidy_hammies %>% 
  glmer(attack_number_of_occurrences ~ sex*treatment+trial_timing+scale(si_cort)
        +weight_at_trial+(1|intruder_id)+(1|juvenile_id),family=poisson,
        glmerControl(optimizer="bobyqa"),data=.)

broom.mixed::tidy(modelA.1) %>% 
  kable("pipe", digits = 3, caption = "**Table 2.** Summary of model A.1, analyzing the interactive effects of maternal treatment and offspring sex on frequency of attack behaviors in the first and last five minutes of the resident-intruder trial.") %>% 
  kable_styling() %>%
  kableExtra::row_spec(.,c(4:6, 10), bold=TRUE) %>% 
  kableExtra::column_spec(column = 1:7, width = "200px")

Table 2. Summary of model A.1, analyzing the interactive effects of maternal treatment and offspring sex on frequency of attack behaviors in the first and last five minutes of the resident-intruder trial.

effect	group	term	estimate	std.error	statistic	p.value
fixed	NA	(Intercept)	0.589	0.706	0.834	0.404
fixed	NA	sexMale	0.430	0.373	1.153	0.249
fixed	NA	treatmentAntibiotic Only	0.604	0.388	1.555	0.120
fixed	NA	treatmentControl	-2.299	1.050	-2.190	0.029
fixed	NA	treatmentStress Only	0.852	0.384	2.221	0.026
fixed	NA	trial_timingLast	-0.321	0.153	-2.103	0.035
fixed	NA	scale(si_cort)	0.023	0.101	0.228	0.819
fixed	NA	weight_at_trial	0.014	0.022	0.645	0.519
fixed	NA	sexMale:treatmentAntibiotic Only	-0.551	0.485	-1.137	0.256
fixed	NA	sexMale:treatmentControl	2.323	1.106	2.101	0.036
fixed	NA	sexMale:treatmentStress Only	-0.784	0.482	-1.626	0.104
ran_pars	juvenile_id	sd__(Intercept)	0.145	NA	NA	NA
ran_pars	intruder_id	sd__(Intercept)	0.000	NA	NA	NA

Model A.1 is a poisson mixed model fitted to predict number of attack bouts.

tidy_hammies %>% 
  ggplot(aes(x=interaction(trial_timing, sex),
             y = attack_number_of_occurrences))+
  geom_boxplot(aes(fill=treatment))+
  labs(x= "Trial timing", 
       y = "Number of attacks", 
       fill = "Maternal treatment")+
  theme(axis.line = element_line(color="black"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank())

Figure 5. Plot of aggression frequency by sex and trial timing.

We then used the lsmeans function from the emmeans R package to conduct pairwise comparisons of the maternal treatment and sex interactions. lsmeans uses a corrected p-value since to minimize the likelihood of a type I error, since we are performing a large number of comparisons. After running the lsmeans, none of our pairwise comparisons are statistically significant using a p-value of 0.05.

#pairwise comparisons with corrected p-value
lsmeans(modelA.1, pairwise~treatment, type = "response") # because there was a significant 2-way interaction, we are more interested in the interaction below.
lsmeans(modelA.1, pairwise~treatment*sex, type = "response") 

lsmeans_output <- lsmeans(modelA.1, pairwise~treatment*sex, type = "response")
lsmeans_contrasts <- lsmeans_output$contrasts
lsmeans_contrasts_df <- as.data.frame(lsmeans_contrasts)
lsmeans_contrasts_df %>% 
  kable("pipe", digits = 3, caption = "**Table 3.** lsmeans comparisons for attack frequency.") %>% 
  kable_styling() %>%
  kableExtra::row_spec(.,14, bold=TRUE) %>% 
  kableExtra::column_spec(column = 1:7, width = "200px")

Table 3. lsmeans comparisons for attack frequency.

contrast	estimate	SE	df	null	z.ratio	p.value
(Antibiotic + Stress Female) / Antibiotic Only Female	0.547	0.212	Inf	1	-1.555	0.777
(Antibiotic + Stress Female) / Control Female	9.964	10.463	Inf	1	2.190	0.358
(Antibiotic + Stress Female) / Stress Only Female	0.427	0.164	Inf	1	-2.221	0.339
(Antibiotic + Stress Female) / (Antibiotic + Stress Male)	0.651	0.243	Inf	1	-1.153	0.945
(Antibiotic + Stress Female) / Antibiotic Only Male	0.617	0.217	Inf	1	-1.370	0.871
(Antibiotic + Stress Female) / Control Male	0.635	0.219	Inf	1	-1.316	0.893
(Antibiotic + Stress Female) / Stress Only Male	0.608	0.231	Inf	1	-1.307	0.897
Antibiotic Only Female / Control Female	18.225	19.168	Inf	1	2.760	0.105
Antibiotic Only Female / Stress Only Female	0.780	0.273	Inf	1	-0.710	0.997
Antibiotic Only Female / (Antibiotic + Stress Male)	1.190	0.445	Inf	1	0.465	1.000
Antibiotic Only Female / Antibiotic Only Male	1.129	0.389	Inf	1	0.352	1.000
Antibiotic Only Female / Control Male	1.161	0.395	Inf	1	0.440	1.000
Antibiotic Only Female / Stress Only Male	1.112	0.418	Inf	1	0.282	1.000
Control Female / Stress Only Female	0.043	0.045	Inf	1	-2.975	0.059
Control Female / (Antibiotic + Stress Male)	0.065	0.067	Inf	1	-2.648	0.139
Control Female / Antibiotic Only Male	0.062	0.064	Inf	1	-2.709	0.120
Control Female / Control Male	0.064	0.066	Inf	1	-2.672	0.131
Control Female / Stress Only Male	0.061	0.063	Inf	1	-2.709	0.120
Stress Only Female / (Antibiotic + Stress Male)	1.525	0.588	Inf	1	1.095	0.958
Stress Only Female / Antibiotic Only Male	1.447	0.485	Inf	1	1.103	0.956
Stress Only Female / Control Male	1.489	0.513	Inf	1	1.155	0.944
Stress Only Female / Stress Only Male	1.425	0.529	Inf	1	0.955	0.980
(Antibiotic + Stress Male) / Antibiotic Only Male	0.949	0.285	Inf	1	-0.176	1.000
(Antibiotic + Stress Male) / Control Male	0.976	0.308	Inf	1	-0.077	1.000
(Antibiotic + Stress Male) / Stress Only Male	0.935	0.299	Inf	1	-0.212	1.000
Antibiotic Only Male / Control Male	1.029	0.307	Inf	1	0.096	1.000
Antibiotic Only Male / Stress Only Male	0.985	0.284	Inf	1	-0.052	1.000
Control Male / Stress Only Male	0.958	0.317	Inf	1	-0.131	1.000

Model assumptions and validation

Below is the code my PI initially used to check model assumptions and check for overdispersion.

# * * Checking model assumptions ------------------------------------------
qqnorm(residuals(modelA.1))
qqline(residuals(modelA.1))
shapiro.test(residuals(modelA.1))


# * * Check overdispersion ------------------------------------------------
overdisp_fun <- function(modelA.1) {
  rdf <- df.residual(modelA.1)
  rp <- residuals(modelA.1,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
overdisp_fun
overdisp_fun(modelA.1)

sum(residuals(modelA.1,type="pearson")^2)/(nrow(tidy_hammies)-length(fixef(modelA.1)-1))

vif(modelA.1)

Using the performance package we learned about in class, I simplified the overdispersion code, and used check_model to visualize vif and overdispersion.

check_model(modelA.1)

Figure 6. Output of performance::check_model function for model A.1.

performance::check_overdispersion(modelA.1)

## # Overdispersion test
## 
##        dispersion ratio =  1.341
##   Pearson's Chi-Squared = 46.932
##                 p-value =  0.086

Aggression Duration

We fitted a linear mixed model for the duration of attack behaviors in the first and last five minutes of the resident-intruder trial. The data were square-root transformed to assume a normal distribution. There were no significant interactions to predict attack behavior duration.

modelA.2 <- lmer(sqrt(attack_total_duration)~ sex+treatment+trial_timing+
                   scale(si_cort)+weight_at_trial+(1|intruder_id)+(1|juvenile_id),
                 data = tidy_hammies)
#summary(modelA.2)
broom.mixed::tidy(modelA.2) %>% 
  kable("pipe", digits = 3) %>% 
  kable_styling() %>%
  kableExtra::column_spec(column = 1:7, width = "200px")

effect	group	term	estimate	std.error	statistic	df	p.value
fixed	NA	(Intercept)	3.608	1.390	2.597	17.398	0.019
fixed	NA	sexMale	0.749	0.452	1.656	17.000	0.116
fixed	NA	treatmentAntibiotic Only	0.324	0.600	0.540	17.000	0.596
fixed	NA	treatmentControl	-0.256	0.586	-0.437	17.000	0.668
fixed	NA	treatmentStress Only	0.397	0.631	0.630	17.000	0.537
fixed	NA	trial_timingLast	-0.009	0.299	-0.031	23.000	0.976
fixed	NA	scale(si_cort)	0.180	0.251	0.717	17.000	0.483
fixed	NA	weight_at_trial	-0.069	0.043	-1.584	17.000	0.132
ran_pars	juvenile_id	sd__(Intercept)	0.677	NA	NA	NA	NA
ran_pars	intruder_id	sd__(Intercept)	0.000	NA	NA	NA	NA
ran_pars	Residual	sd__Observation	1.034	NA	NA	NA	NA

Summary and Conclusions

---

• Overall, both frequency and duration of social behaviors decrease over a 15 minute trial
  
• Maternal treatment has significant effects on attack frequency.
  
• Maternal treatment has sex-specific effects on investigation frequency and duration, escape frequency, and grooming duration.
  
• Trial timing has significant effects on escape duration.
  
• This research suggests that early maternal experience affects behavior during long bouts of stress, and effects can be sex-specific
  
• Mechanisms driving changes in social behaviors during long bouts of stress need to be explored further